#include <stdlib.h>
#include <string.h>

void convert_hetero(char b, char *b1, char *b2) {
    if(b == 'M') {
        *b1 = 'A';
        *b2 = 'C';
    } else if(b == 'Y') {
        *b1 = 'C';
        *b2 = 'T';
    } else if(b == 'S') {
        *b1 = 'C';
        *b2 = 'G';
    } else if(b == 'R') {
        *b1 = 'A';
        *b2 = 'G';
    } else if(b == 'K') {
        *b1 = 'G';
        *b2 = 'T';
    } else if(b == 'W') {
        *b1 = 'A';
        *b2 = 'T';
    } else
        fprintf(stderr, "during covert_hetero %c\n", b);
}
char variant_in_hetero(char b, char ref) {
    char b1, b2;
    convert_hetero(b, &b1, &b2);
    if(b1 != ref && b2 == ref)
        return b1;
    else if(b1 == ref && b2 != ref)
        return b2;
    else {
        fprintf(stderr, "during hetero_var %c\n", b);
        return '\0';
    }
}

#define ALLELE_LEN 500  //allele sequence len
#define SITE_ALLELE_NUM 2000
#define CHR_STR_LEN 100 //for "chr22"

struct sample_site {
    unsigned int chr;//1 2 ... 23(X) 24(Y).
    unsigned int pos;
    unsigned int num_allele;//on this site.
    unsigned char allele_seq[SITE_ALLELE_NUM][ALLELE_LEN];
    unsigned int allele_reads[SITE_ALLELE_NUM];
    unsigned int allele_quality[SITE_ALLELE_NUM];
//quality is sum of base quality in allele in supporting reads, divided by #base in allele after "sam-allele-qual" !.
//for insert, seq is "AT", most than one nt. for delete, seq is "*" and quality from a read is phread_for_delete !!
    struct sample_site *next;
};

//str is "chr22", "chrX", "chrY".
unsigned int get_chr_num(const unsigned char *str) {
    str += 3;
    if(*str == 'X')
        return 23;
    else if(*str == 'Y')
        return 24;
    else
        return atoi(str);
}//OK

int compare_pos(unsigned int chr1, unsigned int pos1, unsigned int chr2, unsigned int pos2) {
    if(chr1 < chr2)
        return -1;
    else if(chr1 == chr2) {
        if(pos1 < pos2)
            return -1;
        else if(pos1 == pos2)
            return 0;
        else
            return 1;
    } else
        return 1;
}//OK

#define ALLELE_INFO_LEN 100000 //"seq:#read:qual/seq:#read:qual/"
unsigned int parse_allele_info(const unsigned char *allele_info, unsigned char allele_seq[][ALLELE_LEN], unsigned int allele_reads[], unsigned int allele_quality[]) {
    unsigned int num = 0;
    for(; *allele_info != '\0'; allele_info = strchr(allele_info, '/') + 1) {
        unsigned char *t = strchr(allele_info, ':');
        int index = t - allele_info;
        strncpy(allele_seq[num], allele_info, index);
        allele_seq[num][index] = '\0';
        t ++;
        allele_reads[num] = atoi(t);
        t = strchr(t, ':') + 1;
        allele_quality[num] = atoi(t);
        num ++;
    }
    return num;
}//OK

//a1,a2, -> array.
unsigned int parse_alleles(const unsigned char *alleles, unsigned char seq[][ALLELE_LEN]) {
    unsigned int index = 0;//as next index.
    unsigned char *comma;
    while(comma = strchr(alleles, ',')) {
        int num = comma - alleles;
        strncpy(seq[index], alleles, num);
        seq[index][num] = '\0';
        alleles = comma + 1;
        index ++;
    }
    return index;
}//OK

unsigned int sum_array(const unsigned int v[], unsigned int num) {
    unsigned int total = 0;
    for(; num > 0; num --)
        total += v[num - 1];
    return total;
}//OK

//file line "chr_num pos allele_info_above"
int to_same_site(FILE *n_file, FILE *t_file, unsigned int *chr, unsigned int *pos, unsigned char n_allele_info[], unsigned char t_allele_info[]) {
    if(feof(n_file) || feof(t_file))
        return 0;

    unsigned int n_chr, n_pos;
    fscanf(n_file, "%u\t%u\t%s\n", &n_chr, &n_pos, n_allele_info);
    unsigned int t_chr, t_pos;
    fscanf(t_file, "%u\t%u\t%s\n", &t_chr, &t_pos, t_allele_info);

    while(1) {
        int i = compare_pos(n_chr, n_pos, t_chr, t_pos);
        if(i > 0) {
            if(feof(t_file))
                return 0;
            else
                fscanf(t_file, "%u\t%u\t%s\n", &t_chr, &t_pos, t_allele_info);
        } else if(i < 0) {
            if(feof(n_file))
                return 0;
            else
                fscanf(n_file, "%u\t%u\t%s\n", &n_chr, &n_pos, n_allele_info);
        } else
            break;
    }
    
    *chr = n_chr;
    *pos = n_pos;
    return 1;
}//OK !!!


//below functions are for test diff significant .........................
void sort_str_value(unsigned char s[][ALLELE_LEN], unsigned int v[], unsigned int n) {
    unsigned int fixed;
    for(fixed = 0; fixed < n; fixed ++) {
        unsigned int moving;
        for(moving = fixed + 1; moving < n; moving ++) {
            if(strcmp(s[moving], s[fixed]) < 0) {
                char t_s[ALLELE_LEN];
                unsigned int t_v;
                strcpy(t_s, s[fixed]), t_v = v[fixed];
                strcpy(s[fixed], s[moving]), v[fixed] = v[moving];
                strcpy(s[moving], t_s), v[moving] = t_v;
            }
        }
    }
}//OK
void move_afterward_from(unsigned char s[][ALLELE_LEN], unsigned int v[], unsigned int n, unsigned int i) {
    for(; n > i; n --) {
        strcpy(s[n], s[n - 1]);
        v[n] = v[n - 1];
    }
}//OK
//on return: correspond allele seq in correspond element, sorted by allele seq, if no such allele, v is set to 0.
unsigned int union_allele(unsigned char s1[][ALLELE_LEN], unsigned int v1[], unsigned int n1, unsigned char s2[][ALLELE_LEN], unsigned int v2[], unsigned int n2) {
    sort_str_value(s1, v1, n1);
    sort_str_value(s2, v2, n2);
    unsigned int i;
    for(i = 0; i < n1 || i < n2; i ++) {
        if(i == n1) {
            strcpy(s1[i], s2[i]);
            v1[i] = 0;
            n1 ++;
        } else if(i == n2) {
            strcpy(s2[i], s1[i]);
            v2[i] = 0;
            n2 ++;
        } else {
            int cmp = strcmp(s1[i], s2[i]);
            if(cmp < 0) {
                move_afterward_from(s2, v2, n2, i);
                strcpy(s2[i], s1[i]);
                v2[i] = 0;
                n2 ++;
            } else if(cmp > 0) {
                move_afterward_from(s1, v1, n1, i);
                strcpy(s1[i], s2[i]);
                v1[i] = 0;
                n1 ++;
            }
        }
    }
    if(n1 == n2)
        return n1;
    else
        return 0;
}//OK !!!

//correspond element in v1 & v2.
double sum_abs_diff_freq(const unsigned int v1[], const unsigned int v2[], unsigned int num) {
    double t1 = sum_array(v1, num);
    double t2 = sum_array(v2, num);
    double stat = 0.0;
    unsigned int i;
    for(i = 0; i < num; i ++) {
        double f1 = v1[i] / t1, f2 = v2[i] / t2;
        if(f1 > f2)
            stat += (f1 - f2);
        else
            stat += (f2 - f1);
    }
    return stat;
}//OK

//sum is fixed.
void random_value(unsigned int v[], unsigned int num) {
    unsigned int total = sum_array(v, num);//remains for rest elements.
    unsigned int i;
    for(i = 0; i < num - 1; i ++) {
        double t = random() / (double) RAND_MAX;
        v[i] = (unsigned int) (t * total);
        total -= v[i];
    }//0 ~ RAND_MAX.
    v[i] = total;
}//OK

